Final Report of NASA Langley Grant NAG- 1-1 145 
Multi-Dimensional High Order Essentially Non-Oscillatory 
Finite Difference Methods in Generalized Coordinates 

Chi-Wang Shu 
Brown University 
June 1, 1990 to May 31, 1998 


/fJ ~ 6 7 


3 


This project is about the development of high order, non-oscillatory type schemes for 
computational fluid dynamics. Algorithm analysis, implementation, and applications are 
performed. Collaborations with NASA scientists have been carried out to ensure that the 
research is relevant to NASA objectives. 

The combination of ENO finite difference method with spectral method in two space 
dimension is considered, jointly with Cai [3]. The resulting scheme behaves nicely for the 
two dimensional test problems with or without shocks. 

Jointly with Cai and Gottlieb, we have also considered one-sided filters for spectral 
approximations to discontinuous functions [2]. We proved theoretically the existence of 
filters to recover spectral accuracy up to the discontinuity. We also constructed such filters 
for practical calculations. 

In [27], jointly with Zang, Erlebacher, Whitaker and Osher, we have applied ENO schemes 
to two and three dimensional compressible Euler and Navier-Stokes equations of gas dynam- 
ics. Practical issues, such as vectorization, efficiency of implementation, cost comparison 
with other numerical methods and accuracy degeneracy effects, are discussed. Numerical 
examples include two and three dimensional compressible homogeneous turbulence, transi- 
tion in a free shear layer, and shock interaction with two dimensional entropy and vorticity 
waves. These cases all have the property that shocks and complicated flow structure co-exist, 
hence non-oscillatory shock capturing and high order accuracy in the smooth part of the flow 
are both important. 

Jointly with Casper and Atkins, we have performed a comparison of the efficiency of the 
two types of ENO schemes. In this paper we have explored the boundary treatment for high 
order ENO schemes in two dimensional problems with a non-rectangular domain. Both wall 
and inflow/outflow boundaries are considered. It is found out that careful treatment of dif- 
ferent boundary conditions is essential for the realization of high order accuracy and stability 
for such problems. We have also found out that the grid orientation effect, which is present 
because of the dimension by dimension feature of the algorithm, is actually diminished with 
the increase of order of accuracy and is, for the test problems we have computed, negligible 
for fourth order schemes. 

Jointly with Perthame, we have studied positivity preserving finite volume schemes in 
one and two space dimensions for arbitrary triangulations. The equations we solve are 
Euler equations of compressible gas, and positivity is preserved for density and pressure. 
Such schemes are weakly stable because, together with conservation, one can prove the L l 
boundedness of the density, momentums, and energy. Such positivity preserving schemes 
also have practical advantages for solving problems near vacuum, for which often negative 
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density and pressure lead to inslability. A general framework and examples are provided in 

[23] 

Among finite difference techniques for flow calculations is the class of compact schemes. 
Compact schemes usually can be designed to have high (spectral like) resolution for various 
waves and are particularly suitable for long time calculation such as direct numerical sim- 
ulation of turbulence. In [5], Cockburn and Shu provided a framework to apply compact 
schemes for shocked or high gradient problems. TVB compact schemes in ID and L <*, stable 
schemes in two and higher dimensions are discussed and tested. 

In [9], E and Shu adapted ENO methodology to solve incompressible flows and investi- 
gated several resolution issues. This is a generalization of the second order Godunov type 
methods used by Bell, Colella and Glaz. The main difference from the compressible flow 
is the presence of the incompressibility condition, which is a global condition and, if not 
treated carefully, may affect local approximation performance. Although the work in [9] 
is still preliminary (only periodic boundary conditions were considered, for which the in- 
compressibility condition can be easily enforced with Fourier transform using the projection 
method), it does provide the key evidence that the high order ENO methodology can be 
applied to incompressible flow and has good resolution power for certain important physical 
quantities. 

The combination of applied analysis and direct numerical simulation to study physics 
of fluids phenomena, is carried out jointly with E. We have studied the effective equations 
and the inverse cascade of energy for the two dimensional Kolmogorov flow in [8], and the 
small scale structure in Boussinesq convection in [10]. Typically, analysis (formal as well 
as rigorous asymptotic analysis) is carried out as far as possible, while carefully designed 
numerical simulation (often with the aid of partial result from analysis) is carried out to 
continue the attack. 

Jointly with Gottlieb, we have been systematically studying the issue of overcoming 
the Gibbs phenomenon, i.e., to recover exponential accuracy from a spectral partial sum 
of a piecewise analytic function. The prototype of recovering in an exponentially accurate 
fashion in the maximum norm, for the Fourier Galerkin case for an analytic but non-periodic 
function, is given in [12]. This paper also sets up the framework for our future developments. 
The resolution issues (number of points per wave) of using Fourier methods for an analytic 
but non-periodic function is studied in [13]. The problem of recovering exponential accuracy 
in a sub-interval, with the assumption that the function is analytic in this sub-interval but 
may have discontinuities at one or both boundary, both for the periodic (Fourier) and for 
the non-periodic (Legendre and Chebyshev) cases, is considered in [14]. The problem of 
recovering exponential accuracy in a sub-interval, with the assumption that the function is 
analytic in this sub-interval but may have discontinuities at one or both boundary, for general 
Gegenbauer polynomial based Galerkin approximations, which includes the Chebyshev and 
Legendre methods as special cases, is investigated in [15]. The problem for the collocation 
case, which is more practical but the proof is much more difficult and involves totally different 
techniques than our proof for the Galerkin case, is contained in [16]. 

Jointly with Wong, we have performed a detailed numerical study about numerical ac- 
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curacy when spectral method is applied to a nonlinear conservation law with discontinuous 
solutions [28]. We assess the accuracy of Fourier Galerkin and collocation method applied to 
Burgers equation with smooth initial condition but with a shock developed in finite time. We 
find that, unlike in the linear PDE case, the moments with respect to analytic functions, in 
particular the first few Fourier coefficients, are no longer very accurate. However the numer- 
ical solution does contain accurate information which can be extracted by a post-processing 
based on Gegenbauer polynomials. 

In [25], which results from an invited talk in the Third International Colloquium of 
Numerical Analysis, we have compared ENO finite difference, finite volume and discontinuous 
Galerkin finite element methods, in terms of their shock resolution, grid orientation effects, 
cost, and ability to handle complicated geometry and boundary conditions. 

Jointly with Harabetian and Osher, we have investigated a novel Eulerian approach for 
simulating vortex motion using a level set regularization procedure [18]. Our approach uses 
a decomposition of the vorticity of the form £ — P(<p)rj , in which both p (the level set 
function) and rj (the vorticity strength vector) are smooth. We derive coupled equations for 
P and 7] which give a regularization of the problem. The regularization is topological and is 
automatically accomplished through the use of numerical schemes whose viscosity shrinks to 
zero with grid size. There is no need for explicit filtering, even when singularities appear in 
the front. The method also has the advantage of automatically allowing topological changes 
such as merging of surfaces. Numerical examples including two and three dimensional vortex 
sheets, two dimensional vortex dipole sheets and point vortices, are given. To our knowledge, 
this is the first three dimensional vortex sheet calculation in which the sheet evolution feeds 
back to the calculation of the fluid velocity. 

Jointly with Jiang, we have been investigating WENO (weighted ENO schemes) [22], 
WENO is a modification and improvement of ENO schemes. Instead of using only one of 
the many candidate stencils based on local smoothness as in ENO, WENO uses a linear 
combination of the contribution from all candidate stencils, each with suitable nonlinear 
weight. The weights are chosen so that in smooth regions, they are close to an optimal 
linear weight which gives the highest possible order of accuracy of an upwind-biased linearly 
stable scheme. Near shocks, however, those stencils which contain the shock are assigned 
weights which are essentially zero. Thus WENO resembles a linear high order upwind biased 
scheme in smooth regions, and ENO near shocks, with a smooth numerical flux function. 
Another advantage of WENO, due to its smoothness of fluxes, is that convergence for Smooth 
solutions can be proven. Also, convergence towards steady states is easier than ENO. 

Jointly with Atkins, we have been performing investigation of the discontinuous Galerkin 
(DG) finite element method for CFD and acoustic calculations in general geometry. In [1], 
we discussed a discontinuous Galerkin formulation that avoids the use of discrete quadrature 
formulas. The application is carried out for one and two dimensional linear and nonlinear 
test problems. This approach requires less computational time and storage than conventional 
implementations but preserves the compactness and robustness inherent in the discontinuous 
Galerkin method. 

Joint work with Erlebacher and Hussaini on shock longitudinal vortex interaction prob- 
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lems, is carried out [11]. In this paper, we have studied the shock/longitudinal vortex in- 
teraction problem in axisymmetric geometry. Linear analysis, shock fitting code, and shock 
capturing ENO are all used, in different parameter regimes, to study various cases of nearly 
linear regime, weakly nonlinear regime, and strong nonlinear regime. Vortex breakdown as 
a function of Mach number ranging from 1.3 to 10 is studied, extending the range of exist- 
ing results. For vortex strengths above a critical value, a triple point forms on the shock, 
leading to a Mach disk. This leads to a strong recirculating region downstream of the shock. 
A secondary shock forms to provide the necessary deceleration so that the fluid velocity can 
adjust to downstream conditions at the shock. 

Jointly with Gottlieb, we have further explored a class of high order TVD (total variation 
diminishing) Runge-Kutta time discretization initialized by Shu and Osher in 1988, which 
is suitable for solving hyperbolic conservation laws with stable spatial discretizations [17]. 
We illustrate with numerical examples that non-TVD but linearly stable Runge-Kutta time 
discretization can generate oscillations even for TVD (total variation diminishing) spatial 
discretization, verifying the claim that TVD Runge-Kutta methods are important for such 
applications. We then explore the issue of optimal TVD Runge-Kutta methods for second, 
third and fourth order, and for low storage Runge-Kutta methods. 

In applications mixed hyperbolic-elliptic systems are sometimes encountered. Examples 
are found in gas dynamics, nonlinear elasticity and Lorenz systems. Numerical solutions of 
mixed type systems often suffer from a lack of guidance in how to stabilize the solution in 
elliptic regions with adequate viscosity. Most of the earlier work in this field use directly 
hyperbolic schemes, e.g., first order Lax-Friedrichs scheme, and second order MarCormack 
scheme. We proposed, in [24], a extension of hyperbolic flux-splitting idea to mixed type 
systems. We have found ways to write a possibly elliptic flux as a sum of two hyperbolic 
fluxes. Our approach works for all 2 x 2 and at least some higher order systems. With this 
flux splitting, any successful hyperbolic numerical operator can be used separately on each 
of the two hyperbolic fluxes. Stability is then obtained as in the hyperbolic case. Numerical 
results in [24] for the van der Waals equation of gas dynamics, using ENO operator for 
each of the split hyperbolic flux, illustrates the good potential for this approach: the phase 
transitions can be correctly resolved and the method can be used as a tool to study the 
evolution of elliptic regions. 

For the discontinuous Galerkin method, the result of Jiang and Shu [20], [21] proves 
a cell entropy inequality for the square entropy, which implies L 2 stability in general and 
convergence in some special cases, for arbitrary triangulations in any space dimension and 
arbitrary order of accuracy, without even using nonlinear limiters. This result is much 
stronger than the corresponding results for finite volume or finite difference schemes, which 
typically require one space dimension, a convex physical flux and no more than second order 
accuracy to prove the same entropy inequality. 

Jointly with Cockburn, we have studied the Local Discontinuous Galerkin methods for 
nonlinear, time-dependent convection-diffusion systems [6]. These methods are an exten- 
sion of the Runge-Kutta Discontinuous Galerkin methods for purely hyperbolic systems to 
convection-diffusion systems and share with those methods their high parallelizability, their 
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high-order formal accuracy, and their easy handling of complicated geometries, for convec- 
tion dominated problems. It is proven that for scalar equations, the Local Discontinuous 
Galerkin methods are L 2 -stable in the nonlinear case. Moreover, in the linear case, it is 
shown that if polynomials of degree k are used, the methods are k - th order accurate for 
general triangulations; although this order of convergence is suboptimal, it is sharp for the 
LDG methods. Preliminary numerical examples displaying the performance of the method 
are shown. 

Again jointly with Cockburn, we wrote the fifth paper in a series [7] in which we con- 
struct and study the so-called Runge-Kutta Discontinuous Galerkin method for numerically 
solving hyperbolic conservation laws. In this paper, we extend the method to multidimen- 
sional nonlinear systems of conservation laws. The algorithms are described and discussed, 
including algorithm formulation and practical implementation issues such as the numerical 
fluxes, quadrature rules, degrees of freedom, and the slope limiters, both in the triangular 
and the rectangular element cases. Numerical experiments for two dimensional Euler equa- 
tions of compressible gas dynamics are presented that show the effect of the (formal) order 
of accuracy and the use of triangles or rectangles, on the quality of the approximation. 

Jointly with Hu, we have presented a discontinuous Galerkin finite element method for 
solving the nonlinear Hamilton-Jacobi equations [19]. This method is based on the Runge- 
Kutta discontinuous Galerkin finite element method for solving conservation laws. The 
method has the flexibility of treating complicated geometry by using arbitrary triangulation, 
can achieve high order accuracy with a local, compact stencil, and are suited for efficient 
parallel implementation. One and two dimensional numerical examples are given to illustrate 
the capability of the method. 

In the lecture note [26], we describe the construction, analysis, and application of ENO 
(Essentially Non-Oscillatory) and WENO (Weighted Essentially Non-Oscillatory) schemes 
for hyperbolic conservation laws and related Hamilton-Jacobi equations. ENO and WENO 
schemes are high order accurate finite difference schemes designed for problems with piece- 
wise smooth solutions containing discontinuities. The key idea lies at the approximation level, 
where a nonlinear adaptive procedure is used to automatically choose the locally smoothest 
stencil, hence avoiding crossing discontinuities in the interpolation procedure as much as 
possible. ENO and WENO schemes have been quite successful in applications, especially 
for problems containing both shocks and complicated smooth solution structures, such as 
compressible turbulence simulations and aeroacoustics. This lecture note is basically self- 
contained. It is oui hope that with this note and with the help of the quoted references, the 
reader can understand the algorithms and code them up for applications. Sample codes are 
also available from the author. 

About high older essentially non-oscillatory (ENO) finite difference schemes, jointly with 
Zeng, we have applied ENO method to the viscoelastic model with fading memory [29] The 
memory term is treated by introducing new variables and rewrite the system by adding 
more differential equations but without explicit memory terms. The appearance of the 
memory terms regularizes the solution somewhat, and in many cases it is still a theoretically 
open question whether shocks will develop from smooth initial data. We have performed 



theoretical analysis about the linearized system for large time, and have applied ENO scheme 
to study the nonlinear system for both local time and large time. The high order accuracy 
and sharp, non-oscillatory shock transition allow us to obtain fine resolution for tens of 
thousands of time steps, and to study the shock interactions after the formation of shocks. 

There are 29 papers published in this period acknowledging support from this NASA 
grant. These papers are listed in the References. 
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